home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Monster Media 1996 #15
/
Monster Media Number 15 (Monster Media)(July 1996).ISO
/
math
/
alged34.zip
/
ALGEDSRC.ZIP
/
ALGEXPR.C
< prev
next >
Wrap
C/C++ Source or Header
|
1996-06-06
|
20KB
|
787 lines
/*--------------------------------------------------------------------
Alged: Algebra Editor henckel@vnet.ibm.com
Copyright (c) 1994 John Henckel
Permission to use, copy, modify, distribute and sell this software
and its documentation for any purpose is hereby granted without fee,
provided that the above copyright notice appear in all copies.
*/
#include "alged.h"
/*--------------------------------------------------------------------
associate
rotate the association on add,sub,mul,div
a+b-c => a-c+b => b+a-c => b-c+a
note, we make use of the fact that ADD+1 = SUB and ADD is even.
*/
void associate(node *p) {
node *a,*b;
int opr;
opr = p->kind;
b = p->rt;
if (opr==DIV && b->kind==MUL) { /* special handle for bisected */
p->rt = b->rt;
b->rt = b->lf;
b->lf = p->lf;
p->lf = b;
b->kind = DIV;
b = p->rt;
}
if (opr==ADD || opr==MUL || opr==SUB || opr==DIV) {
a = p;
while ((a->lf->kind|1) == (opr|1)) {
a->kind = a->lf->kind;
a->rt = a->lf->rt;
a = a->lf;
}
a->kind = opr;
if (opr&1) { /* subtr or divide */
a->rt = b;
}
else {
a->rt = a->lf;
a->lf = b;
}
}
else if (opr==EQU) { /* commute equality */
p->rt = p->lf;
p->lf = b;
}
}
/*--------------------------------------------------------------------
commute
commute on add,sub,mul,div
e.g. a/b ==> b^(-1)/a^(-1)
*/
void commuteNOTUSED(node *p) {
node *a,*b;
a = p->lf;
b = p->rt;
if (p->kind==MUL || p->kind==ADD || p->kind==EQU) {
p->lf=b; p->rt=a;
}
else if (p->kind==DIV || p->kind==SUB) { /* non-commutative */
if (a->kind==p->kind+1 &&
b->kind==p->kind+1 &&
a->rt->kind==NUM &&
b->rt->kind==NUM) { /* has coefficients */
p->lf=b; p->rt=a;
a->rt->value = -a->rt->value;
b->rt->value = -b->rt->value;
}
else { /* need to make coefficients */
p->lf = newoper(p->kind+1);
p->rt = newoper(p->kind+1);
p->lf->lf = b;
p->lf->rt = newnum(-1);
p->rt->lf = a;
p->rt->rt = newnum(-1);
}
}
}
/*--------------------------------------------------------------------
a^x * a^y ==> a^(x+y)
x^a * y^a ==> (xy)^a
*/
int expjoin(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=expjoin(p->parm[i]);
a = p->lf;
b = p->rt;
if (p->kind==MUL || p->kind==DIV) {
if (a->kind==EXP && b->kind==EXP &&
equal(a->lf,b->lf)) { /* a^x * a^y = a^(x+y) */
freetree(a->lf);
a = a->rt;
freenode(p->lf);
b->kind = p->kind-2;
p->kind = EXP;
p->lf = b->lf;
b->lf = a;
}
else if (b->kind==EXP &&
equal(a,b->lf)) { /* a * a^y = a^(1+y) */
freetree(a);
a = newnum(1);
b->kind = p->kind-2;
p->kind = EXP;
p->lf = b->lf;
b->lf = a;
}
else if (a->kind==EXP &&
equal(a->lf,b)) { /* a^x * a = a^(x+1) */
freetree(b);
b = newnum(1);
a->kind = p->kind-2;
p->kind = EXP;
p->lf = a->lf;
a->lf = a->rt;
a->rt = b;
p->rt = a;
}
else if (a->kind==EXP && b->kind==EXP &&
equal(a->rt,b->rt)) { /* x^a * y^a = (xy)^a */
freetree(b->rt);
p->rt = a->rt;
a->rt = b->lf;
freenode(b);
a->kind = p->kind;
p->kind = EXP;
}
else return r;
++r;
}
return r;
}
/*--------------------------------------------------------------------
remove sub and div
x/y = x*y^-1
x-y = x+y*-1
*/
int nosubdiv(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=nosubdiv(p->parm[i]);
a = p->lf;
b = p->rt;
if (p->kind==SUB || p->kind==DIV) {
++r;
--p->kind;
if (b->kind==NUM && p->kind==ADD) /* a-2 = a+(-2) */
b->value = -b->value;
else if (b->kind==NUM && !whole(b->value))
b->value = 1.0/b->value;
else {
a = newoper(p->kind+2); /* MUL or EXP */
a->lf = b;
a->rt = newnum(-1);
p->rt = a;
}
}
return r;
}
/*--------------------------------------------------------------------
bisect node
This converts expressions to canonical form in which
1. + * are left assoc, ^ is right assoc.
2. scope of / is increased, e.g. (a/b)*c ==> ac/b,
3. x+y*-2 is changed to x-y*2
3' x+(-2) is changed to x-2
4. x*y^-2 is changed to x/y^2
4' x*(0.5) is changed to x/2
5. scope of ^ is reduced
The name "bisect" means that each MUL clag is broken into two pieces the
numerator elements and the denominator ones.
*/
int bisect(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=bisect(p->parm[i]);
a = p->lf;
b = p->rt;
switch (p->kind) {
case ADD:
case MUL:
/*--------------------------------------------------------------------
Add or Multiply
*/
if (b->kind==p->kind || /* a+(b+c) = a+b+c */
b->kind==p->kind+1) { /* a+(b-c) = a+b-c */
swingb;
i = b->kind;
b->kind = p->kind;
p->kind = i;
}
else if (p->kind==ADD && /* a+(-2) = a-2 */
b->kind==NUM &&
b->value < 0) {
p->kind = SUB;
b->value = -b->value;
}
else if (b->kind==p->kind+2 && /* a+b*-2 = a-b*2 */
b->rt->kind==NUM &&
b->rt->value < 0) {
++p->kind;
b->rt->value = -b->rt->value;
}
/*--------------------------------------------------------------------
Add only
*/
else if (p->kind==ADD &&
aop(b->kind)) { /* a+(b-c) = (a+b)-c */
swingb;
p->kind = b->kind;
b->kind = ADD;
}
/*--------------------------------------------------------------------
Multiply only
*/
else if (p->kind==MUL &&
a->kind==p->kind+1) { /* a-b+c = a+c-b */
p->rt = a->rt;
a->rt = b;
++p->kind;
--a->kind;
}
else if (p->kind==MUL &&
a->kind==p->kind+2 && /* b*-2+a = a-b*2 */
a->rt->kind==NUM &&
a->rt->value < 0) {
++p->kind;
p->lf = b;
p->rt = a;
a->rt->value = -a->rt->value;
}
else break;
++r; break;
case SUB:
case DIV:
/*--------------------------------------------------------------------
Subtract or divide
*/
if (b->kind==p->kind+1 && /* a-b*-2 = a+b*2 NOT NECES. */
b->rt->kind==NUM &&
b->rt->value < 0) {
--p->kind;
b->rt->value = -b->rt->value;
}
/*--------------------------------------------------------------------
Subtract only
*/
else if (p->kind==SUB && /* a-(-2) = a+2 */
b->kind==NUM &&
b->value < 0) {
p->kind = ADD;
b->value = -b->value;
}
else if (p->kind==SUB &&
aop(b->kind)) { /* a-(b+c) = (a-b)-c */
swingb;
p->kind = (ADD+SUB) - b->kind;
b->kind = SUB;
}
/*--------------------------------------------------------------------
Divide only
*/
else if (p->kind==DIV &&
b->kind==p->kind) { /* a-(b-c) = a+c-b */
p->lf = b;
p->rt = b->lf;
b->lf = a;
--b->kind;
}
else if (p->kind==DIV &&
a->kind==p->kind) { /* a-b-c = a-(b+c) */
swinga;
--a->kind;
}
else break;
++r; break;
case EXP:
/*--------------------------------------------------------------------
exponent
*/
if (a->kind==EXP) { /* (x^y)^z = x^(y*z) */
swinga;
a->kind = MUL;
}
else if (b->kind==NUM && /* a^(-2) = 1/a^2 */
b->value < 0) {
p->kind = DIV;
p->lf = newnum(1);
p->rt = newoper(EXP);
p->rt->lf = a;
p->rt->rt = b;
b->value = -b->value;
}
else break;
++r; break;
}
return r;
}
/*--------------------------------------------------------------------
exponent expand - expand any integer exponents less than 100.
*/
int exexpand(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=exexpand(p->parm[i]);
a = p->lf;
b = p->rt;
if (p->kind==EXP &&
b->kind==NUM && /* a^n = a^(n-1)*a */
b->value > 1 &&
b->value <= maxpow &&
whole(b->value)) {
if (--b->value == 1) {
p->lf = deepcopy(a);
freenode(b);
}
else {
b = newoper(EXP);
b->lf = deepcopy(a);
b->rt = p->rt;
p->lf = b;
}
p->rt = a;
p->kind = MUL;
++r;
}
return r;
}
/*-----------------------------------------------------------------
within - if p is a factor within q, then return the rest of q.
if the result is not null, then it is a new copy.
*/
node *within(node *p, node *q) {
node *a=q,*b=NULL;
while (a->kind==MUL) {
if (equal(p,a->rt)) {
if (b) {
b->lf = a->lf; // skip over match
q = deepcopy(q); // make copy
b->lf = a; // repair
}
else q = deepcopy(a->lf);
return q;
}
else if (equal(p,a->lf)) {
if (b) {
b->lf = a->rt; // skip over match
q = deepcopy(q); // make copy
b->lf = a; // repair
}
else q = deepcopy(a->rt);
return q;
}
b = a;
a = b->lf;
}
return NULL;
}
/*--------------------------------------------------------------------
ComDeno - find common denominators.
*/
int comdeno(node *p) {
int i,r=0;
double x;
node *a,*b,*nu,*de;
for (i=0; i<p->nump; ++i)
r+=comdeno(p->parm[i]);
a = p->lf;
b = p->rt;
if (aop(p->kind) && (a->kind==DIV || b->kind==DIV)) {
if (a->kind==DIV && b->kind==DIV) {
if (equal(a->rt,b->rt)) { /* a/x + b/x = (a+b)/x */
freetree(b->rt);
de = a->rt;
a = a->lf;
b = b->lf;
freenode(p->lf);
freenode(p->rt);
}
else if (!!(nu=within(a->rt,b->rt))) { /* a/x + b/xy = (ay+b)/xy */
freetree(a->rt);
a = cons(nu,MUL,a->lf);
de = b->rt;
b = b->lf;
freenode(p->lf);
freenode(p->rt);
}
else if (!!(nu=within(b->rt,a->rt))) { /* a/xy + b/x = (a+by)/xy */
freetree(b->rt);
b = cons(nu,MUL,b->lf);
de = a->rt;
a = a->lf;
freenode(p->lf);
freenode(p->rt);
}
else { /* a/c + b/d = (ad+bc)/cd */
a->kind = b->kind = MUL;
de = b->rt;
b->rt = a->rt;
a->rt = de;
de = newoper(MUL);
de->lf = deepcopy(b->rt);
de->rt = deepcopy(a->rt);
}
}
else if (a->kind==DIV) { /* a/b + c = (a + bc)/b */
a->kind = MUL;
de = a->lf;
a->lf = b;
b = a;
a = de;
de = deepcopy(b->rt);
}
else { /* a + b/c = (ac + b)/c */
b->kind = MUL;
de = b->lf;
b->lf = a;
a = b;
b = de;
de = deepcopy(a->rt);
}
nu = newoper(p->kind);
nu->lf = a;
nu->rt = b;
p->kind = DIV;
p->lf = nu;
p->rt = de;
++r;
}
return r;
}
/*-----------------------------------------------------------------
distribute2
This applies the distributive laws
1. division over addition.
*/
int distribute2(node *p) {
node *a,*b,*c;
int i,r=0;
for (i=0; i<p->nump; ++i)
r+=distribute2(p->parm[i]);
i=0;
c = p->parm[i];
if (p->kind==DIV && aop(c->kind)) {
b = deepcopy(p->parm[1-i]);
a = newoper(p->kind);
a->parm[i] = c->parm[1-i];
a->parm[1-i] = b;
c->parm[1-i] = p->parm[1-i];
p->parm[1-i] = a;
p->kind = c->kind;
c->kind = a->kind;
++r;
}
return r;
}
/* This is used to control the direction of distribute */
static int topdown = 0;
/*-----------------------------------------------------------------
distribute
This applies the distributive laws
1. multiplication over addition.
2. (ab)^2 = a^2*b^2
3. x^(a+b) = x^a*x^b
*/
int distribute(node *p) {
node *a,*b,*c;
int i,r=0;
if (!topdown)
for (i=0; i<p->nump; ++i)
r+=distribute(p->parm[i]);
for (i=0; i<2; ++i) {
c = p->parm[i];
if (p->kind==MUL && aop(c->kind) ||
i && p->kind==EXP && aop(c->kind) ||
!i && p->kind==EXP && (c->kind==MUL || c->kind==DIV)) {
b = deepcopy(p->parm[1-i]);
a = newoper(p->kind);
a->parm[i] = c->parm[1-i];
a->parm[1-i] = b;
c->parm[1-i] = p->parm[1-i];
p->parm[1-i] = a;
p->kind = c->kind;
c->kind = a->kind;
if (i && a->kind==EXP) p->kind+=2; /* change ADD to MUL */
++r;
break; /* just to be safe */
}
}
if (topdown && !r)
for (i=0; i<p->nump; ++i)
r+=distribute(p->parm[i]);
return r;
}
/*-----------------------------------------------------------------
distribute_c
this applies the distributive law over multiplication which is
not governed by another multiplication or equality.
*/
int distribute_c(node *p) {
int r=0;
topdown = 1;
if (p->kind==EQU) {
r+=distribute_c(p->lf);
r+=distribute_c(p->rt);
}
else if (p->kind==MUL) {
r+=distribute(p->rt);
r+=distribute_c(p->lf);
}
else r+=distribute(p);
topdown = 0;
return r;
}
/*--------------------------------------------------------------------
fixassoc
+ * are left assoc, ^ is right assoc.
*/
int fixassoc(node *p) {
int i,r=0;
node *a,*b;
for (i=0; i<p->nump; ++i)
r+=fixassoc(p->parm[i]);
a = p->lf;
b = p->rt;
switch (p->kind) {
case ADD:
case MUL:
if (b->kind==p->kind || /* a+(b+c) = a+b+c */
b->kind==p->kind+1) { /* a+(b-c) = a+b-c */
swingb;
i=b->kind; b->kind=p->kind; p->kind=i;
}
else break;
++r; break;
case SUB:
case DIV:
if (b->kind==p->kind) { /* a-(b-c) = a-b+c */
swingb;
--p->kind;
}
else if (b->kind==p->kind-1) { /* a-(b+c) = a-b-c */
swingb;
++b->kind;
}
else break;
++r; break;
case EXP:
if (a->kind==EXP) { /* (x^y)^z = x^(y*z) */
swinga;
a->kind = MUL;
}
else break;
++r; break;
}
return r;
}
/*-----------------------------------------------------------------
get term, return p without coefficient
*/
node *get_term(node *p, int oper, double *r) {
*r = 1.0; /* default coeff */
if (p->kind==oper && p->rt->kind==NUM) {
*r = p->rt->value;
return p->lf;
}
return p;
}
/*--------------------------------------------------------------------
Combine all terms with common base within a clag.
*/
int combine(node *p) {
int i,oper,r=0;
node *a,*b;
node *t1,*t2;
double c1,c2;
for (i=0; i<p->nump; ++i)
r+=combine(p->parm[i]);
a = p->lf;
b = p->rt;
oper = p->kind;
if (oper!=ADD && oper!=MUL) /* not a clag */
return r;
i = 1;
if (a->kind!=oper) {
a=p; i=0;
}
t1 = get_term(a->parm[i],oper+2,&c1);
t2 = get_term(b,oper+2,&c2);
if (equal(t1,t2)) { /* same base, combine terms */
if (a->parm[i]==t1)
freetree(t1);
else { /* free one of the expressions */
freetree(b);
b = a->parm[i];
t2 = t1;
}
c2 += c1;
if (b==t2) { /* make room for coeff */
b = newoper(oper+2);
b->lf = t2;
b->rt = newnum(0);
}
b->rt->value = c2;
if (c2==1.0) { /* remove unit coeff */
freenode(b->rt);
freenode(b);
b = t2;
}
if (i) { /* cleanup... 2 stage */
nodecpy(p,a);
freenode(a);
p->rt = b;
}
else { /* 1 stage */
nodecpy(p,b);
freenode(b);
}
++r;
}
return r;
}
/*-----------------------------------------------------------------
Simplify
This function converts a rational expression to normal form.
Some of the attributes of normal form: no negative exponents,
terms are sorted over mul and add.
The first movenums pushes all numbers to the left so that
when the sortnode runs it pushes them to the right and combines
them.
*/
void simplify(node *p) { simplify2(p,0); }
void simplify2(node *p, int slow)
{
while (calcnode(p,1)) if (slow) slow=debug(p);
while (fixassoc(p)) if (slow) slow=debug(p);
while (nosubdiv(p)) if (slow) slow=debug(p);
while (fixassoc(p)) if (slow) slow=debug(p);
while (movenums(p,0,MUL)) if (slow) slow=debug(p);
if (sortnode(p,MUL) && slow) slow=debug(p);
while (movenums(p,0,ADD)) if (slow) slow=debug(p);
if (sortnode(p,ADD) && slow) slow=debug(p);
while (combine(p)) if (slow) slow=debug(p);
while (calcnode(p,1)) if (slow) slow=debug(p);
if (sortnode(p,MUL) && slow) slow=debug(p);
if (sortnode(p,ADD) && slow) slow=debug(p);
while (bisect(p)) if (slow) slow=debug(p);
while (calcnode(p,1)) if (slow) slow=debug(p);
while (movenums(p,0,MUL)) if (slow) slow=debug(p);
}
/*-----------------------------------------------------------------
substitute -- prereq tgt MUST be an EQU.
*/
void substitution(node *p) {
int i;
if (equal(p,tgt->lf))
movenode(p,deepcopy(tgt->rt));
else for (i=0; i<p->nump; ++i)
substitution(p->parm[i]);
}
/*-----------------------------------------------------------------
insert key
has special handling for equations...
*/
void insertkey(int opr) {
node *tmp,*t1,*t2;
if (!src || !tgt) return;
if (tgt->kind==EQU) {
t1 = tgt->lf; t2 = tgt->rt;
} else t1 = t2 = tgt;
if (src->kind==EQU) {
tmp = src->lf;
src->lf = newoper(opr);
src->lf->lf = tmp;
src->lf->rt = deepcopy(t1);
tmp = src->rt;
src->rt = newoper(opr);
src->rt->lf = tmp;
src->rt->rt = deepcopy(t2);
src->rt->lf = tmp;
}
else if (opr==EXP) {
tmp = newnode();
nodecpy(tmp,src);
src->kind = EXP;
src->nump = 2;
src->lf = cons(tmp,EXP,deepcopy(t1));
src->rt = cons(newnum(1),DIV,deepcopy(t2));
src = src->lf;
}
else {
tmp = newnode();
nodecpy(tmp,src);
src->kind = opr ^ 1;
src->nump = 2;
src->lf = newoper(opr);
src->rt = deepcopy(t2);
src = src->lf;
src->lf = tmp;
src->rt = deepcopy(t1);
}
}
/*-----------------------------------------------------------------
invert left side of equality, i=0 or 1
*/
void cross_eq(node *p, int i) {
node *t;
if (p->kind != EQU) return;
t = p->lf;
if (t->kind==ADD ||
t->kind==SUB ||
t->kind==MUL ||
t->kind==DIV) {
p->lf = t->parm[i];
t->parm[i] = p->rt;
p->rt = t;
if (i && !(t->kind&1)) { /* swap lf and rt */
t = t->rt;
p->rt->rt = p->rt->lf;
p->rt->lf = t;
t = p->rt;
}
if (!i || !(t->kind&1)) { /* toggle operator */
t->kind ^= 1;
}
}
}